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,  ABSTRACT 

//■  -  ■  ' 

Considers^  direct  iterative  method  for  solving  the  linear  system  AU  =  Y 
which  arises  from  the  discretization  of  a  boundary  value  problem  involving  an 
elliptic  partial  differential  operator  L  of  order  2m.  Such  iterative 
methods  split  A  into  a  difference  A  =  M  -  N.  The  critical  question  is 
the  behavior  of  p  =  max{|\|,  X  an  eigenvalue  of  In  this  report  we 

develop  a  theory  for  the  determination  of  the  asymptotic  behavior  of  p  as 

D 

h  -+■  0.  This  theory  depends  on  recognizing  (h  N)  as  a  weak  approximation 
to  a  differential  operator  q  of  order  _<  m.  Then,  under  "appropriate  con- 

p 

ditions",  p  ~  1  -  AQh  where  AQ  is  the  real  part  of  the  minimal  eigen¬ 
value  of  Lu  =  Xqu .  The  study  of  the  appropriate  conditions  leads  one  to 
three  interesting  mathematical  questions:  (i)  relevance,  (if)  Spectral 
Approximation  and  (iii)  general  estimates. 
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1 .  Introduction 

An  important  problem  that  arises  often  in  applications  is  the 
solution  of  an  elliptic  boundary-value  problem 


(1.1) 

Lu  =  f  in  ft 

(1.2) 

Bu  =  0  on  3ft 

where  ft 

operator 

simplest 

is  a  bounded  (or  unbounded)  smooth  domain  in  IR^  .  The  boundary 
B  is  in  general  a  linear  differential  operator.  However,  the 

case  is  (for  second  order  problems)  merely 

(1.2*) 

u  =  0  on  3ft  . 

Such  elliptic  boundary-value  problems  are  frequently  stated  in  their 

"weak  form."  Let  H  (ft)  be  the  functions  defined  on  ft  which  have  m'th 
m 

derivatives  in  L-(ft).  Let  H  c  H  (ft)  be  a  subspace  that  incorporates 
i.  mm 

(in  an  appropriate  manner)  the  boundary  conditions  (1.2).  Let  a(u,<M  be 

a  bilinear,  continuous  and  coercive  form  defined  on  H  ®H  : 

m  m 

(1.3) 

a(u,4>)  is  bilinear,  so  that 

(1.3a) 

a(au1  +  6u  2 ,  d> )  =  aa(u1t$)  +  &a(u2»<)>)  and 

(1.3b) 

a(u,a<j>1  +S<J>2)  =  aa(u,4>-|)  +  6a (u ,4>2 )  ; 

(1.4) 

a(u,<J>)  is  continuous,  so 

there  is 

a  constant  K  >  0  such  that 

(1.4a) 

1  a (u ,4>)  |  1  K ||  u  || m  •  ||  4>|lm  »  Vu,  4>  e  Hm  ; 

L 
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at 


(1.5)  a(u,<(>)  is  coercive,  i.e., 

a(u,u)  is  real  and  there  is  a  constant  KQ  >  0  such  that 
(1.5a)  Kq i |  u  ||^  <  a(u,u),  Vu  e  Hm  . 

Then,  under  certain  appropriate  conditions,  the  problem  (1 .1 )-(l .2)  takes 
the  form: 


pind  u  e  H  so  that 
-  m - 

(1.6)  a(u,<|>)  =  F($) ,  ¥0  e  Hm 

where  F  is  a  linear  functional  defined  on  H  . 

m 

Condition  (1.4)-(1.6)  guarantees  that  there  is  a  unique  solution  u*. 
However,  in  general,  one  cannot  explicitly  exhibit  u*.  Then,  the  desire 
for  precise  quantitative  information  about  u*  leads  tc  numerical  approaches. 
The  problem  is  "discretized,"  i.e.,  the  problem  is  set  in  a  finite-dimensional 
subspace,  and  we  obtain  a  system  of  linear  algebraic  equations  of  the  form 

(1.7)  AU  =  F  . 

Here  A  is  an  nxn  matrix,  n  is  large,  and  (somehow) 

(1.8)  A  ~  a  . 

One  approach  to  the  solution  of  (1.7)  is  a  direct  iterative  method  based 
on  a  "splitting"  of  the  matrix  A.  We  write 

(1.9)  A  =  M  -  N  , 


whence  (1.7)  takes  the  form 
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(1.10)  MU  =  NU  +  F  . 

Given  a  first  guess,  say  U°,  we  then  obtain  iterates  UJ  from  the 
scheme 

(1.11)  MUj+1  =  NUj  +  F  . 

We  are  now  almost  ready  to  discuss  the  problem  of  this  report  -  the 
asymptotic  behavior  of  the  spectral  radius  of  the  i teration  matri x 

(1.12)  K  =  M_1N  . 

It  is  well  known  that  the  convergence  of  the  iterative  scheme  (11)  is 
determined  by  the  spectral  radius  of  K  .  Therefore  we  are  concerned 
with  the  eigenvalue  problem 

(1.12a)  AMX  =  NX,  X  t  0 

and  especially  with 

(1.12b)  p:=  max  | X J  . 

We  suppose  that  this  "splitting"  has  been  arranged  in  some  regular  way 
and  OUR  PROBLEM  is  the  study  of  the  ASYMPTOTIC  BEHAVIOR  OF  p  as 

n  ■*  <*>  . 

This  is  an  old  problem.  There  are  results  of  Shortley  and  Weller  -  1938 
[18],  a  beautiful  theory  by  Young  -  1951  [21],  which  extended  work  of  Frankel 
1950  [7],  an  intriguing  paper  by  Garabedian  -  1956  [8],  a  fine  book  by  Varga  - 
1962  [19];  and  right  now  there  is  an  active  school  in  Brussels,  Belgium 
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working  with  Beauwens  [2],  [3],  [4]  applying  methods  of  analysis  to 
these  problems.  Of  course,  many,  many  others  (at  least  a  zillion)  whom 
I  have  not  mentioned  have  made,  and  are  making  contributions.  I  apologize 
to  all  whose  names  have  been  omitted. 

The  theory  of  Young  is  based  on  a  combinatorial  concept  -  Property  A. 
The  reader  need  not  be  expert  in  this  theory.  However,  I  do  wish  to 
mention  one  interesting  and  important  consequence  of  this  theory. 

For  (block)  Jacobi  schemes  which  have  Property  A,  if  A  is  an 
eigenvalue  of  M-1N,  so  is  (-X)! 

Our  approach  to  this  problem  can  be  outlined  as  follows.  We  rewrite 
(1.12a)  as 

A(M-N)X  =  (l-X)NX , 

and  if  A  |  0  (A  =0  is  a  particularly  uninteresting  case)  then 

0.13)  AX  =  [ijlA][hPN]X  . 

Ahp 

That  is, 

(1.14a)  AX  =  pNX  , 

(1.14b)  p  =  ~~  ,  N  =  hpN, 

Ahp 

where  p  is  as  yet  undetermined.  Because  A  arises  from  L  in  some 
natural  way  we  hope  to  find  a  related  eigenvalue  problem  of  the  form 


1. 
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(1.15)  LU  =  Aqu  ,  Bu  =  0  on  3ft  , 

where  q  is  a  differential  operator  of  degree  lower  than  2m  (<  m  in  fact!). 
Thus  in  some  sense 

(1.16)  hpN  ~  q  . 


Then,  "if  all  goes  right,"  there  will  be  an  eigenvalue  X  of  (1.12a) 

with 

|A|  =  p 

and 


1-X 

Ahp 


=  y  = 


+  o(l). 


where  Am  =  AQ  +  IT  is  a  particular  eigenvalue  of  (1.15).  Then 


X  =  — —  =  1  -  A  hp  +  o(hp) 

1 +uhp  m 

and 

p  =  1  -  Aghp  +  o(hp)  . 

In  his  basic  paper  Garabedian  [8]  began  with  equation  (1.11) 

k+1  k+1 

as  applied  to  the  error  e  -  u  -  u  and  wrote 


M( 


j 

-) 


This  expression  led  him  to  a  time-dependent  problem 

Jyt 


Y  ll  =  auxt  +  Suvt  +  Lu  ’ 


which  in  turn  led  to  the  eigenvalue  problem 


(1.18)  L<J>  =  A[yu  -au  "  8u  ]  ,  Bu  =  0 

x  y 

and  the  condition 

(1.19)  o  =  e  U  =  1  -  AQAt  +  Q(A td) 

where  AQ  is  the  "smallest"  eigenvalue  of  (1.8). 

As  one  develops  this  idea  one  finds  many  interesting  mathematical 
questions  that  must  be  resolved  before  one  can  assert  the  validity  of 

(1.19) .  Garabedian  neither  raised  nor  resolved  these  questions.  Yet  in 
the  particular  cases  he  discussed,  we  believe  his  results  are  correct. 

In  one  case  (the  five  point  star  for  the  Laplace  operator)  these  results 
follow  independently  from  the  theory  of  Young  [21],  [22].  There  is  a 
large  amount  of  computational  experience  that  indicates  the  validity  of 

the  results  in  the  other  case  (the  nine-point  star  for  the  Laplace  operator). 

In  this  report  we  develop  the  theory  in  a  very  general  setting. 

The  theory  and  results  are  useful  and  interesting  in  themselves.  However, 
we  shall  also  emphasize  three  interesting  mathematical  questions  that 
this  theory  brings  to  the  forefront.  These  are: 

(i)  Spectral  approximation. 

(ii)  Relevance:=  a  study  of  properties  of  those  eigenvalues  A 
wi  th 

|  A  l  =  p  . 

(iii)  Estimates:3  clearing  out  the  unimportant  and  providing  a 
basis  for  proofs  of  Spectral  approximation. 

I 

i 


i  *=TT  I  J 
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2.  Examples 


We  start  with  a  simple  1-dimensional  example 


(2.1) 

(2.2) 


-u"  =  f  ,  0  <  x  <  1 

u(0)  =  u(l)  =  0  . 


Let  us  consider  the  simplest  finite  difference  approach  to  this  problem. 


(2.3) 


and  let  A  be  the  n*n  tridiagonal  matrix 


2  -1 
-1  2 


O 


(2.4) 


A  =  h 


=  £  L-1,2,-1] 


o 


Most  people  would  expect  to  find  (-L)  rather  than  -jj-  .  However,  we 

n  n 

have  chosen  this  normalization  both  because  it  is  appropriate  for  the 
presentation  of  our  theory  and  because  this  normalization  is  appropriate 
from  the  finite  element  point  of  view. 


V 
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To  this  end,  with  each  vector 

U  —  (u-j  ,^2*  —  un) 

we  associate  a  continuous  piecewise  linear  function 

u(x,h) 

which  is  linear  on  each  interval  kh  <  x  £  (k+l)h  and  satisfies 

u(kh,h)  =  u^  .  u(0,h)  =  u(l,h)  =  0  . 

The  space  of  all  such  functions  is  designated  .  Then  (see  Figure  1) 
for  every  two  such  vectors  U,  V  we  have 


(2.5) 


V*AU 


1 

v'  (x.h)u' (x,h)dx  =  a(u,v)  . 

0 


Thus  in  this  example  we  see  exactly  how  A  and  a(u,4>)  are  associated 
recall  (1.8). 


Figure  1 
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In  any  case,  the  boundary  value  problem  (2.1)-  (2.2 
the  linear  system 

(2.6)  AU  =  F 

where  the  vector  F  is  related  to  the  function  f(x)  . 
in  the  simplest  case 

F  =  h(f(h),f(2h),...f(nh))T  . 

The  (point)  Jacobi  iterative  scheme  (which  we  certainly 
but  discuss  as  a  simple  example)  is  given  by 


(2.7) 


0  1 
1  0 


0 

0 


1 

0 


0  0 
1  0 
0  1  0 


1 


0 


Hence 


with 

(2.8a) 


is  replaced  by 


For  example, 


do  not  recommend. 
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and 


It  happens  that  the  eigenvalues  and  eigenvectors  of  are 

known  [19].  These  are 

(2.9)  =  cosirkh  ,  Xk  =  (simrkjh) . 


n 


Remark  2:  Consider  the  eigenvalues  X  with 

|X|  =  p  . 

There  are  two,  X ,  and  X  .  In  this  case 
1  n 

(2.10)  p  =  X.|  =  cos  7rh  =  1  -  iT^h^  +  o(h^)  . 

Let  us  write 

(2.11 )  A=|l-W-WT  =  D-  W-WT 

where  W  is  a  strictly  lower  triangular  matrix.  Then  the  Jacobi  scheme 
just  discussed  is  given  by 

(2.12)  |  I  Uk+1  =  (W+WT)Uk  +  F 
and  the  corresponding  eigenvalue  problem 

XMX  =  NX  ,  X  f  0 

is  a  "nice"  standard  problem  for  the  eigenvalues  of  a  real  symmetric  matrix 
N  relative  to  a  positive  definite  matrix  M. 

The  associated  point  Gauss-Seidel  iteration  scheme  is  given  by 

(2.13)  (|  I-  W)Uk+1  =  WTUk  +  F 

and  its  associated  eigenvalue  problem  is  not  a  regular  problem:  N  is  not 
symmetric  and  M  is  not  positive  definite.  But  there  are  some  saving 
features.  The  matrix 


K  =  M_1N 
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is  a  nonnegative  irreducible  matrix.  Hence  the  Perron-Frobenius  theory 
[19]  tells  us  that  there  is  exactly  one  eigenvalue  A  with 

m  -  p 


and,  in  fact,  that  A  satisfies 


A  =  p  ! 

Indeed,  the  theory  of  Young  [2l],  [22]  tells  us  that 

(2.14)  p  =  (cos  Trh ) ^  =  1  -  ir^h^  +  O(h^)  . 

Another  important  scheme  based  on  the  representation  (2.11)  is  the 
SOR-w  (successive  overrelaxation)  scheme.  Let  u>  be  chosen  with 

0  <  u  <  2  . 


Consider  the  scheme 

(2.15)  l(D-wW)Uk+1  =  l[wWT  +  (l-u>)D]Uk  +  F  . 

This  method  is  known  to  be  convergent  [19]  and  if  y  f  0  is  an  eigen¬ 
value  of  this  SOR-w  scheme  then  there  is  a  k,  1  _<  k  £  n,  for  which  y 

the  equation 

(2.16)  (y  +w  -1)^  =  u^(cos  nkh)^y  . 

For  our  purposes  we  need  only  observe  that  (see  page  111  of  [19]):  for 

each  uj,  0  <  u>  <  2,  the  solutions  y  are  either  real  and  positive  or 

complex,  and  both  y  and  y  are  eigenvalues.  In  either  case,  there  is 


satisfies 
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an  eigenvalue  X  with 


M  =  P 


and 

(2.17)  Re  X  >  0  . 

Remark:  In  all  three  cases  we  have  discussed  the  "relevance"  question  - 
as  defined  in  section  1.  In  a  few  pages  we  shall  see  why  this  is  indeed 
relevant. 

Let  us  now  return  to  our  general  problem 
(*)  XMU  =  NU  . 

Then 


X(M-N)U  =  XAU  =  (l-X)NU  . 
Hence  if  X  ^  0  we  obtain 


AU 


[hpN]U  . 


That  is 

(2.18a)  AU  =  yNU 

where 

(2.18b)  u  =  —  ,  N  =  hpN 

Xhp 
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The  exponent  p  has  not  yet  been  chosen.  Consider  the  Jacobi 
iterative  scheme  (2.7). 

Let  us  write  (2.18a)  in  the  "weak  form":  Find  U  so  that 
(2.18c)  V*AU  =  uV*NU  ,  VV  e  Cn  . 


Using  the  basic  (finite  element)  fact  (2.5)  we  have 
1 


(2.19) 


0 


D  1  * 

v'(x,h)u'(x,h)dx  =  yhP  . 


Since  the  left-hand-side  of  (2.19)  is  an  integral  and  the  right-hand-side 
is  "almost"  a  Riemann  Sum,  we  choose 


(2.20) 


P  =  2  . 


We  write 


Ub+l”2u.+u.  i 

twvj  ■  2“k  *  ht~  -h-  1 


and  (2.19)  becomes 


(1  +  h2p)  V*AU  =  2y  h£  vfcu  , 


or 


(1  +h2u) 


1  1 
v' (Xjh)u' (x,h)dx  ~  2p  j  v(x,h)u(x,h)dx  . 
0  0 


2 

If  we  ignore  the  term  h  p  we  see  that  we  have  convinced  ourselves 
that  the  eigenvalue  M  of  (2.18)  -  with  p  =  2  -  is  related  to  eigen¬ 
values  of  the  weak  problem 
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(2.21) 


1 

a(u,4>)  =  A 

0 


2u(x)<{>(x)dx  . 


Our  first  question  is:  does  (2.21)  have  eigenvalues?  Because  (2.21) 
is  the  weak  form  of 

(2.22a)  -u"  =  Aq(x)u  ,  0  <  x  <  1 

(2.22b)  u(0)  =  u(l)  =  0 

wi  th 

(2.22c)  q(x)  =  2 

the  answer  is  yes.  In  particular,  we  come  to  the  first  part  of  our 
discussion  of  spectral  problems.  The  problem  (2.21)  (or  equivalently 

(2.22) )  possesses  a  minimal  eigenvalue  AQ  >  0  .  In  fact 

(2.23)  A0=^tt2. 

Our  next  question  is:  Do  there  exist  eigenvalues  of  (2.18)  near  Ag? 
Let  us  be  more  precise:  Let  us  define 

Spectral  Approximation  (a):  Let  <5  >  be  given.  Then  there  is  an  hg  >  0 
such  that  for  each  h,  0  <  h  _<  hg,  there  is  an  eigenvalue  u(h)  of  (2.18) 
satisfying 


I Ag-u(h) |  <  6  . 

Now  our  question  is:  Is  Spectral  Approximation  (a)  true?  In  this  case  the 
answer  is  yes! 
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Because  the  answer  is  yes,  there  is  an  eigenvalue 


X  =  A(h)  =  - 7 


of  (*)  and 


1+h  y(h) 


1  ?  ?  ? 

X  =  1  -  j  TT  il  +  o(h^)  . 


Hence,  by  the  definition  of  p 
(2.24) 


p  _>  1  -  ir^h2  +  o(h2)  . 


In  order  to  complete  the  discussion  and  reobtain  (2.10)  we  must  know  that 
the  eigenvalue  problems  (2.18)  and  (2.21)  are  relevant  to  the  eigenvalue 
problem  (*).  For  example,  suppose  we  did  not  know  the  eigenvalues  of  (2.18) 
and  (as  is  not  the  case)  the  eigenvalues  X  with  j X [  =  p  all  satisfied 


X  <  0  ; 


then  the  corresponding  p  satisfy 

|u(h)  I  =  l^rl  f  +  00  • 

All 

Thus,  these  eigenvalues  of  (2.18)  are  not  relevant  to  the  limiting  eigen¬ 
value  problem  (2.21).  Let  us  define 

Relevance:  There  are  constants  Cg  >  0  ,  hg  >  0,  and  for  0  <  h  <  hg 
there  is  an  eigenvalue  X  with  | X J  =  p  and 


VJ 


(2.25) 
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In  addition  to  the  question  of  relevance  we  must  discuss 

Spectral  Approximation  (b):  Let  p(hn)  be  a  sequence  of  eigenvalues 
of  (2.18)  with  hn  -*■  0  .  Suppose 

U  (  hn  )  -  u  . 

n 

Then  p^  is  an  eigenvalue  of  the  limiting  equation  (2.21). 

In  this  case,  because  of  the  variational  characterization  of  eigenvalues 
of  (2.18)  and  the  variational  characterization  of  the  eigenvalues  of  (2.21), 
it  is  not  difficult  to  prove  that  Spectral  Approximation  (b)  is  true  and  we 
can  complete  the  argument  as  follows.  Let  X ( h )  be  a  relevant  eigenvalue 
with 


A(h)  =  p(h)  . 


Because  we  know  (from  very  general  considerations)  that  the  Jacobi  scheme  is 
convergent,  we  deduce  from  (2.24)  that 

1  -  \  Tr2h2  +  o ( h2 )  <  x(h)  <  1  , 

and  hence  the  relevance  condition  is  satisfied  with 


C 


0 


Therefore  the  associated  eigenvalues  p(h)  of  (2.18)  are  uniformly  bounded. 
Hence  a  subsequence,  say  p(h  ),  converges  to  an  eigenvalue  of  (2.21),  say  p^ 
Then 


^  =  •*.*0(1) 
h 
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and 

A  =  1  -  u  h2  +  0(h2)  . 

Hence,  because  u  is  of  the  form 

00 

-  1  vK1 

^co  '  2  k  * 

we  get 

2 

(2.26)  A  =  p  <  1  --y7r2h2  +  o(h2)  . 

Combining  (2.26)  with  (2.24)  yields  (2.10). 

This  approach  was  used  by  Parter  [11],  [12]  to  establish  a  general  theorem  - 
in  the  self-adjoint  case.  More  recently  Parter  and  Steuerwalt  [14]  have  extended 
that  argument  to  finite  difference  equations  for  the  general  second  order  ellip¬ 
tic  operator,  and  to  related  parabolic  problems. 

Now  let  us  consider  the  Gauss-Seidel  iterative  scheme.  Following  our  recipe 
we  obtain  (2.18)  with 


(2.27)  N  =  l  [0,0,1] 
and 

(2.28)  V*NU  =  hP_1  l  vkuk+1  . 

k-1 


JljJL 


Hence,  arguing  as  before  we  have 


(2.29a)  P  =  2 

and 

(2.29b)  V*NU  =  h  l  v  u  +  h[h  £  v  (  ^  k)]  . 

k=l  k  h  k=l  K  n 

Because 

our  equation  (2.19)  becomes 
1 

(2.30)  |  v'(x)u'(x)dx  =  Mh^vkuk  +  uhe(v.u) 

0 

where  (using  standard  inequalities) 

(2.30b)  |e(v,u)J  <  CJ|v  ||  -)|  u'  ||  .  . 

u  l2  L2 

In  this  case  we  are  led  to  the  related  problem 

1 

(2.31)  a(u,<j>)  =  A  j  u(x)i(x)dx  ,  V4>  eS^  . 

0 

This  is  the  weak  form  of  the  problem 
(2.32a)  -u"  =  Aq(x)u  ,  0  <  x  <  1 

(2.32b)  u(0)  =  u(l)  =  0 

with 

(2.32c)  q(x )  =  I  . 


•i - - 
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In  this  case 


and  we  are  led  to  conjecture  that 

(2.33)  p  =  1  -  u2h2  +  o(h2)  . 

As  we  know  -  see  (2.14)  -  this  is  correct!!  Let  us  see  what  we  must  prove 
in  order  to  apply  our  general  arguments. 

Once  more,  general  considerations  imply  that  the  method  is  convergent. 
And,  as  we  have  indicated  earlier,  the  Perron-Frobeni us  theory  implies  that 
p  is  itself  an  eigenvalue.  However,  as  we  have  said,  in  this  instance  the 
eigenvalue  problem  (1.12a)  is  not  a  self-adjoint  problem  while  the  eigenvalue 
problem  (2.32a)  is  a  self-adjoint  problem.  Hence  in  this  case  the  verifica¬ 
tion  of  the  spectral  approximation  property  is  not  so  simple. 

Now  let  us  consider  the  point  SOR  iterative  scheme  (2.15).  With  our 
knowledge  of  the  nature  of  the  final  results  based  on  the  SOR  theory  of 
Young  [21],  [22]  let  us  set 

(2.34a)  u  =  2  -  Ch  , 

(2.34b)  N  =  hN  ,  p  =  ^  , 

i.e., 

P  =  1  • 


V 


A  direct  computation  shows  that 
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(2.35) 


V*NU  =  h£[^C]vkuk  +  h[vk[-k^-  -—I 
+  0(h2Ivkuk)  . 


Thus  we  are  led  to  compare  the  basic  discrete  eigenvalue  problem 
(2.18a)  with  the  eigenvalue  problem 


(2.36a)  -u“  =  A[-^  Cu+u']  ,  0  <  x  <  1 

(2.36b)  u(0)  =  u(l)  =  0  . 


Thus,  in  this  case 

(2.37)  1=TC’lr- 

Once  more  we  ask  our  basic  questions: 

(i)  Does  (2.36)  have  eigenvalues? 

(ii)  Is  Spectral  Approximation  (a)  true? 

(iii)  Is  Spectral  Approximation  (b)  true? 

(iv)  What  about  the  relevance  of  the  eigenvalues  of 
the  discrete  problem  (2.18)  to  the  eigenvalue 
problem  (2.36)? 

In  this  simple  case  we  know  all  the  answers  and  may  apply  our  theory. 

We  hope  these  examples  have  given  some  notion  of  how  one  proceeds  with 
these  ideas.  However,  in  all  the  examples  so  far,  (because  we  have  been 
dealing  with  "point"  iterative  schemes)  it  would  seem  that  the  "weak  form" 
of  eigenvalue  problem  is  unnecessary.  We  can  readily  recognize  (2.18a)  as 
a  consistent  approximation  to  the  eigenvalue  problem 
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(2.38)  -u"  =  Aqu  ,  u(0)  =  u(l  )  =  0  . 

Hence  we  close  this  section  with  a  simple  -  but  complex  -  example:  the 
"4-point  Jacobi  Scheme."  However,  let  us  once  more  remind  the  reader  that 
this  procedure  is  not  being  suggested  as  a  means  of  solving  this  simple  prob¬ 
lem.  It  is  an  interesting  example  of  the  ideas  involved  in  this  theory.  Let 

(2.39)  n  =  4a 
and  let  A^  be  the  4x4  matrix 

(2.40a)  A4  -l  C-1,2,-1]  , 

and  give  M  and  N  the  induced  block  tridiagonal  structures 


(2.40b) 


A 


4 


0  0  0  0 
0  0  0  0 

where  R  = 

0  0  0  0 
.10  0  0 

As  usual,  our  iterative  scheme  is  given  by  (1.11).  We  take 
(2.41a)  N  =  h2N  ,  p  =  2  , 


(2.41b) 

Then 
(2.42) 
Clearly,  if 


v4s+l  *  v4s  *  u4s+l  *  u4s 

then  (2.42)  suggests  that  V*NU  is  a  Riemann  sum  based  on  a  Ax  of  4h 
But  the  multiplier  of  (2.42)  is  h  .  So  we  rewrite  (2.42)  as 


(2.43a) 


V*NU  *  4h  l  v4sU4Sq 
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where 

(2.43b)  q  =  \  =  f  . 

Thus  we  are  led  to  the  eigenvalue  problem 
(2.44a)  -u"  =  Aqu  ,  u(0)  =  u(l )  =  0  , 

(2.44b)  q  =  \  . 

In  this  way,  after  answering  our  mathematical  questions  of  Spectral 
Approximation  and  "relevance,"  we  have 

(2.45)  p4  =  1  -  2tt2h2  +  o(h2)  . 

Similarly,  if  we  consider  the  k-point  block  Jacobi  iterative  scheme  we 
obtain 

(2.46)  P|<  =  1  -  |  ttV  +  o(h2)  . 

While  the  results  (2.45)  and  (2.46)  are  interesting,  we  invite  the 
reader  to  see  N  -  as  given  by  (2.40c)  *  as  a  consistent  approximation 
to  the  function  q(x)  =  1/2  .  In  fact,  the  operators  N  that  arise 
in  block  schemes  with  large  blocks  all  behave  as  does  (2.40c).  They 
are  rather  nasty  operators.  They  are  the  "homogenization"  operators  of 
Babuska  [l ]  and  Lions  [5]. 
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3.  A  More  Formal  Outline 

Let  us  reconsider  what  we  have  done.  Beginning  with  a  matrix  A 
that  is  related  to  the  differential  operator  L  or,  equivalently, 
to  the  bilinear  form  a ( u , <f> )  ,  we  write 

A  =  M  -  N  ,  (M""'  exists ) 

and  consider  the  eigenvalue  problem 

(3.1)  AMU  =  NU  . 

Our  concern  is  with  the  spectral  radius 

(3.2)  p  =  maxlA]  . 

We  rewri te  (3.1 )  as 

(3.3a)  V*AU  =  yV*NU  ,  ¥V  e Cn  , 

where 

(3.3b)  y  =  —  ,  N  =  hPN  . 

Ahp 

Having  done  this  we  seek  to  discover  an  operator  q  such  that 

q:  Hm  L - 

and,  in  some  way,  we  may  recognize  (3.3)  as  an  approximation  to  the 
weak  eigenvalue  problem 


(3.4) 


a(u,<J>)  =  A  $q[u]dxdy  ,  V<J>  c  H  . 

J  J  ••• 
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Then  we  assume 

A .1 )  The  eigenvalue  problem  (3.4)  possesses  a  "minimal " 


eigenvalue 


(3.5) 


A  =  An  +  iT  . 
m  u 


with  AQ  >  0  .  By  "minimal"  we  mean  that  for  any 
eigenvalue  A  we  have 


0  <  Aq  <  Re  A  , 
and  if  Re  A  =  AQ  then 


|A„I  i  1 A | 


A. 2)  Spectral  Approximation  (a),  (b):  for  h  small  enough 
there  is  an  eigenvalue  u(h)  near  A^  ,  and  the  limit 
of  every  convergent  sequence  of  eigenvalues  { n ( hh ) } 
of  (3.3a)  is  an  eigenvalue  A  of  (3.4). 


A. 2)  Relevance;  There  are  positive  constants  Cg  ,  hg  so 
that  for  every  h  <  hg  there  is  an  eigenvalue  \  of 
the  basic  eigenvalue  problem  (3.1)  that  satisfies 

|A]  =  p 


-  c° ' 

Under  these  circumstances  it  is  an  easy  matter  to  show  that 


p  =  1  -  Aghp  +  o(hp) 


(3.6a) 


When  one  can  prove  the  truth  of  Spectral  Approximation  (a),  (b)  but 
cannot  discuss  the  relevance  question  one  obtains 


(3.6b)  p  >  1  -  AQhp  +  o(hp)  . 

This  result  involves  two  of  the  three  mathematical  topics  discussed 
in  the  introduction.  We  have  not  yet  discussed  “estimates." 

In  the  discussion  of  the  point  Jacobi  iterative  scheme  (2.7)  we  were 
led  to  (2.19).  We  then  wrote 


(3.7a)  h l  vk(u(<+1+uk_1  )  =  h£(2)vkuk  +  e(v,u) 

and 


1 


(3.7b)  hi  vk(uk+1+uk-1)  =  j 

0 

In  dealing  with  the  4-point  Jacobi 


2v(x)u(x)dx  +  E(v,u)  . 

scheme  we  wrote  -  see  (2.42) 


a-1 

(3.8a) 

hJ1(''4stlU4s*54sU4stl) 

=  hi  2v4su4s+  e(v,u) 

and 

a-1 

1 

=  |  (^)vudx  +  E(v,u) 
0 

(3.8b) 

hs=/V4s+lU4s  +V4sU4s+1  ^ 

In  these  formulae,  the  terms  e(v,u),  2(v,u),  E(v,u),  E(v,u)  represent 
"errors"  that  can  be  estimated.  The  question  is,  is  it  sufficient 
merely  to  recognize  these  errors  as  being  small  provided  that  u(x),  v(x) 
are  "nice,  smooth"  functions,  or  are  more  precise  estimates  required?  In 
work  to  date  we  have  met  two  situations.  In  the  first  case  -  and  by  far 


our 
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the  most  common  one  —  one  can  easily  prove:  if  (3.3a)  holds  and 
v(x)  is  smooth  then 

(3.9a)  V*AU  =  a(u,v)  +  o(  ||  v  ||s  +||  u )]  H  ) 

m 

(3.9b)  V*NU  =  [  vq[u]dxdy  +  o(  ||  v  ||  +||  u  ||  H  ) 

where  ||  v||  denotes  an  appropriate  norm  of  v  .  In  this  case  the 
proof  of  the  validity  of  the  Spectral  approximation  hypothesis  is 
relatively  straightforward.  In  the  finite  element  case  this  result 
follows  from  the  (by  now)  standard  theories  -  e.g.  see  [10].  For  the 
finite  difference  case  we  carried  out  an  appropriate  modification  of 
the  standard  argument  for  a  special  problem  in  the  non  self-adjoint 
case  [13].  The  ideas  in  that  work  generalize.  Of  course,  when  both  (3.3a) 
and  (3.4)  are  self-adjoint  the  Spectral  Approximation  problem  is  much  simpler 
One  can  use  the  variational  principles  that  characterize  eigenvalues  to  ob¬ 
tain  the  desired  results.  Such  proofs  of  Spectral  Approximation  based  on 
variational  arguments  are  quite  old  -  see  [20]. 

In  our  general  work  with  Steuerwalt  on  finite  element  equations  we 
used  very  precise  estimates,  which  we  will  discuss  in  Section  5,  to  obtain 
a  new  convergence  theorem.  Moreover,  as  we  will  see  in  Section  5,  those 
estimates  are  extremely  useful  in  the  discovery  of  "q." 


I 
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4.  Finite  Difference  Methods  -  2nd  Order  Problems 

We  now  consider  a  special  case  that  arises  often  in  practice:  2nd 
order  linear  elliptic  operators  in  2  dimensional  regions,  i.e. 


(4.1) 


Lu:=  '  ^(aux}x  +(bux)y  +^buy^x  +{cUy^y^ 


+  d^u 


x 


+  d2uy  +  dQu  . 


Actually,  the  fact  that  we  are  considering  only  2-dimensional  problems  is 
merely  a  notational  convenience.  The  ideas  of  this  section  generalize 
to  elliptic  problems  in  any  number  of  dimensions. 

In  many  finite  difference  equations  for  such  problems  we  may  write 
the  unknown  as  ukj. ,  where  we  imagine 

ukj  “  “(kAx.jAy)  • 


That  is,  the  underlying  space  of  functions  are  continuous  piecewise 
linear  functions  and  the  unknowns  {u.  .}  are  merely  the  function  values 

K  J 

at  selected  points  in  the  region  Q.  In  this  setting  it  is  convenient 
to  use  double  indices  to  describe  the  vector  (u^.)  of  unknowns. 
Correspondingly,  the  matrix  A  of  our  basic  problem  (7)  is  described  by 
pairs  of  double  indices.  We  write  (7)  as 


(4.2) 


op 


a,  .  u 
kj,ap  ap 


Let's  assume  that 


V*AU  =  a(u,v)  +  e  (u,v) 

3 


T 


(4.3a) 
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where 

(4.3b)  ea(u ,v)  =  o(a(u,u)  +a(v,v))  . 

In  many  iterative  schemes  we  have  a  relatively  simple  splitting  (1.9) 
that  satisfies 


Condition  S: 

(1)  If 

(2)  If 


a.  . 
kj.crp 

n.  . 
kj, ay 


=  0 
t  0 


then 

then 


nkj, ay 
nkj ,ay 


m.  .  =  0  . 

kj,  ay 


-a 


kj  ,ay 


If  Condition  S  is  satisfied  then  we  take 


(4.4)  N  =  h2N  ,  y  -  ^—4  . 

Xlv 

Remark:  The  fact  that  p  =  2  and  we  use  h  is  related  to  the  fact 
that  this  is  a  second  order  equation,  not  to  the  fact  that  this  problem 
is  set  in  2  dimensions.  In  d  dimensions  we  have  (we  hope  the  notation 
is  clear) 


ak  k  k  a  a  a  =  °^d"2)  • 

f  *2 *  *  *  a  *1  2*  *  *  Q 


Let  us  compute  V*NU  .  We  have 


(4.5) 

Thus 


V*NU  =  h Z  i  l  n.  . 


V., 


kj  ay 


kj ,ou  kj  ay  * 


(4.6a)  V*NU  -  h  Ul  nkj>0y]vkjukj  *  ^(u,v) 


where 
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-1 


(4.6b)  J(u.v)  -  b2  I  "kj,ou(%-“kj)J  • 

In  finite  difference  methods  there  is  usually  a  fixed  integer  r  >  0  such  that 


unless 


nkj, ay 


=  0 


|  k-a |  <  r  and  |j-y|  <  r  . 


Hence  for  "smooth"  functions  u  we  expect  e(u,v)  to  be  small. 
Specifically,  we  anticipate 

(4.6c)  |e(u,v)|  _<  Kh{h2[|vf(j|2}15  ^(u.u)*5  . 

Thus,  at  first  glance,  our  candidate  for  q  is  given  by 


I  " 

ay 


kj.au 


However,  as  we  have  seen  in  the  case  of  the  4-point  block  Jacobi  scheme, 
it  may  well  be  that  qRj  =  0  for  most  (k,j)  .  Indeed,  if  our  "blocks" 
are  large,  this  will  always  be  the  case.  Hence,  we  must  take  this  fact 
into  account. 

Assuming  that  the  nonzero  values  of  qkj-  occur  with  some  regular 
pattern  throughout  the  mesh  we  may  try 


"(vV  ~“5m 

where  a  is  a  parameter  that  describes  the  density  of  points  (k,j) 
at  which  qkj  is  (structurally)  nonzero. 
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Now  one  has  a  candidate  for  q  and  must  proceed  to  consider  the 
problems  of  Spectral  approximation  and  Relevance.  As  mentioned  earlier, 
the  work  of  Parter  and  Steuerwalt  [14]  and  Parter  [12]  discuss  this 
problem  in  great  detail. 


1 


4 
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5.  Finite  Element  Equations 

In  dealing  with  finite  element  equations  one  can  be  overwhelmed  with 
the  complexity  of  the  search  for  "q."  In  a  recent  paper  [15]  Parter  and 
Steuerwalt  presented  a  general  theory  for  the  linear  systems  that  arise 
in  the  finite  element  approximation  of  elliptic  problems. 

In  addition  to  presenting  this  general  approach  in  the  finite  element 
setting,  [15]  deals  with  a  specific  group  of  problems.  These  problems 
illustrate  the  complexity  as  well  as  the  value  and  significance  of 
appropriate  "estimates."  Let 

(5.1 )  Q  :=  { (x,y),  0  <  x,y  < 1 } 

and  let  L  be  given  by  (4.1).  We  assume  this  problem  is  coercive.  The 
finite  element  space  is  the  space  of  tensor  products  of  Hermite  cubic 
splines  based  on 

1  .  1 

4x  "  fTT  •  4y  *  5“fT  ■ 

Mx  Ky 

Thus  to  each  point  (kAx.jAy)  we  associate  a  4-vector 

ukj  '  Cukj-Ax(ux)kj-<4y(“y)kr‘Sx<',<uxy)y]  • 

Therefore  in  general  at  each  point  there  are  144  coefficients  and  V*NU  involves 

v,  v  ,  v  ,  v  ,  u,  u  ,  u  ,  u  .  In  Figure  3  we  exhibit  the  "9-point  block  star" 

*  y  *  y  *y 

associated  with  the  simplest  problem 

-Au  =  f  ,  Ax  =  Ay  =  h  . 

Since  the  finite  element  theory  usually  gives  estimates  on  the  norm 


A. 


Figure  3„  The  nine-point  decomposition  of 
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and  does  not  estimate  the  point  values  of  u  ,  u  and  u  it  would  seem 

x  y  xy 

difficult  to  discover  the  needle  "q"  in  that  messy  haystack.  However, 
the  finite  element  spaces  generally  have  some  regularity  and  estimates. 

In  this  case  we  obtain 

Theorem  5.1 :  Let  h  =  /Ax Ay  .  Then  there  is  a  constant  K  >  C,  depending 
on  r  =  ,  such  that  for  every  u  ,  v  e  Sh  and  every  <J>  e  C^(fi)  we  have 

(5.2a)  h2  I|v  .|2  <  K||v||f  , 

i  j  J  2 

(5.2b)  h2  l  u.jV.^.j  =  J  J  uvpdxdy  +  <5(u,v,<J>)  , 

^  J 

(5.2c)  h2  I  i  Kh2t|Vv|!^  . 

(5. 2d)  h2  l  l  (Ax)  1(Ay)CX2|(Dotu)  ,i2  <  Kh2||Vu||f  , 

ij  lotl-1  10  “  L2 

' 

(5.2e)  fr  [  £  (Ax)  Z(A  y)  Z|  (D°u) .  .(DBv).  .  |i<  Kn(u,v,h)  , 

ij  Ial+IBI>1  J  1J  “ 

.  j 

where 

(5.3a)  |6(u,v,<f>)|  <_  K[1  +  ||  V^H^jnfu+v.u-v.h) 

and 

(5.3b)  n(u,v,h)  =  h[  ||  u  || ,  ||  Vv  || ,  +  |jv||  ||Vu||  +  h||Vu||  ||Vu||  ]  . 

l2  2  2  2  2  2 
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Let  us  see  what  these  estimates  mean.  The  estimate  (5.2b)  asserts 

that 


h2i  ;ijuu  '  1 1 5udxdy  • 

The  estimate  (5.2c)  justifies  the  sort  of  index  shifting  that  we  did  in 

(4.5),  (4.6).  Generally  speaking,  (5.2c),  (5. 2d),  (5.2e)  tell  us  that  we 

may  ignore  many  terms  -  both  in  the  discovery  of  q  and  in  verification 

* 

of  Spectral  approximation.  This  is  particularly  true  if  the  splitting 
satisfies  Condition  S  of  section  4. 

Moreover,  with  the  estimates  of  Theorem  5.1  one  can  prove:  There  is  a 
relevant  eigenvalue  A  with  | X J  =  p  that  satisfies  (A. 3)  provided  that 
there  is  such  an  eigenvalue  with  eigenvector  U  satisfying 

Re  U*NU  >  0  . 

See  Theorem  5.1  of  [15].  Notice  that  block  Jacobi  schemes  that  possess 
Property  A  always  possess  this  property.  This  result  yields  new  convergence 
proofs  that  do  not  depend  on  having  self-adjoint  problems  or  equations  of 
"positive  type."  The  main  results  are: 

Theorem  5.2:  Consider  the  k-line  block  Jacobi  method  (horizontal  lines). 
Then 

(5.4)  q(x,y)  =  -y-  c(x,y)  ,  r  =  ^  , 

and  the  spectral  radius  p,(k)  satisfies 


L  _ 
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(5.5a)  Pj(k)  =  1  -  ||  rQ(Ay)2  +  o(h2) 

where  Tg  is  the  minimal  eigenvalue  of  the  eigenvalue  problem 

(5.5b)  Li^  =  A  c(x,y)iji  . 

The  spectral  radius  Pg,.(lO  Tor  the  k-line  Gauss-Seidel  scheme 
satisfies 

(5.6a)  Pg<j(k)  =  1  -  |  krQ(Ay)2  +  0(h2)  . 

Similarly,  the  spectral  radius  p  (k)  of  the  SOR  scheme  with  optimal 
choice  of  u>  satisfies 

(5.6b)  Pb(k)  -  1  -  2(^  Tg^Ay  +  o(h)  . 

Theorem  5.3.  Let  the  unknowns  be  ordered  lexicographically  and  consider 
the  point  Gauss-Seidel  iterative  scheme.  Then 

q(*,y)  =  |t|  Oa(x,y)  +~  c(x,y)]  , 

and 

P  >  1  -  AQh2  +  o(h2)  . 

Notice ;  In  this  case  we  cannot  claim  asymptotic  equality.  In  the  self- 
adjoint  case  we  know  [19»  Theorem  3.6]  that  the  method  is  convergent.  After 
all,  Gauss-Seidel  is  always  convergent  for  positive  definite  problems.  Unfor¬ 
tunately,  it  appears  that  none  of  the  proofs  of  this  fact  gives  any  quali¬ 
tative  information  about  the  eigenvalues  X  with  |X|  =  p  . 
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6.  Finite  Elements  and  Point  SOR 

The  results  of  the  previous  section  do  not  yield  results  for  the  "point" 
SOR  scheme  applied  to  the  finite  element  equations.  Nevertheless  there  are 
several  reasons  for  seeking  such  results.  For  one  thing,  experimental  re¬ 
sults  by  Fix  and  Larsen  [6]  and  Rice  [17]  imply  that  for 

h  <  T5 

the  SOR  iterative  schemes  are  preferable  to  direct  solution  methods. 

Part  of  the  strength  of  the  finite  element  method  is  the  fact  that  one 
need  not  have  a  nice  regular  array  of  mesh  points.  One  puts  points  where 
they  are  needed  and  effects  the  necessary  triangulation.  Thus  we  would 
not  always  expect  to  be  able  to  use  nice  regular  blocks  for  iterative  methods. 
Hence  in  these  situations  we  would  expect  to  be  using  point  SOR  rather  than 
k-line  schemes  or  other  regular  block  schemes. 

Finally,  because  this  point  SOR  scheme  does  not  fall  within  the  Young 
theory  (i.e.,  we  do  not  have  property  A)  there  does  not  exist  an  efficient 
adaptive  algorithm  (see  [9])  to  determine  the  optimal  to. 

For  these  reasons  Parter  and  Steuerwalt  undertook  a  study  of  point  SOR 
for  the  model  problem 


-Au  =  f 

where  the  coefficients  of  A  are  given  in  Figure  3.  We  set 

N  =  hN  ,  to  =  2  -Ch 


and  find  that 
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(6.4) 
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c=v 

and  Aq(c)  is  linear  in  c  for  0  <  c  <  . 

While  we  cannot  assert  that  we  have  found  p(u>),  we  have  found  a  candi¬ 
date  for  a  lower  bound  p(w)  that  exhibits  the  characteristic  features  of  p(w) 
in  the  classical  case  where  property  A  holds.  It  is  our  hope  that  the 
nature  of  p(w)  will  be  useful  in  the  construction  of  adaptive  algorithms 
for  the  determination  of  an  optimal  w. 


■k  V‘> 
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7„  Concluding  Remarks 

In  the  previous  sections  we  have  attempted  to  elucidate  a  theory 
for  the  estimation  of  the  asymptotic  behavior  of  the  spectral  radii 
of  iterative  schemes  for  elliptic  problems.  This  theory  has  two  main 
interesting  mathematical  points. 

(1)  The  theory  does  not  lose  sight  of  the  elliptic  operator  L. 

(2)  The  theory  leads  one  to  study  other  mathematical  questions 
that  are  interesting. 

Finally,  while  we  have  not  attempted  to  describe  the  usefulness 
of  the  theory,  it  is  very  useful  -  even  though  the  basic  constant  AQ 
is  generally  not  accessible. 
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